## Agroforestry Paper Analysis import
## tree dataset
library(readxl)
nyando1 <- read_excel("C:/Users/LORERO/OneDrive - CGIAR/Desktop/ICRAF/Current Projects/Nyando Project/Data/2019 Nyando Data.xlsx", 
    sheet = "set")
View(nyando1)

nyando1$Block = as.factor(nyando1$Block)
nyando1$type_r = as.factor(nyando1$type_r)
nyando1$agr = as.factor(nyando1$`Agroforestry Practice`)

# Create subset for Total Lower
# Project
nyandolp <- subset(nyando1, type_r == 
    "Lower Project")
View(nyandolp)

# Create subset for Total Lower
# Control
nyandolc <- subset(nyando1, type_r == 
    "Lower Control")
View(nyandolc)

# Create subset for Total Middle
# Project
nyandomp <- subset(nyando1, type_r == 
    "Middle Project")
View(nyandomp)

# Create subset for Total Middle
# Control
nyandomc <- subset(nyando1, type_r == 
    "Middle Control")
View(nyandomc)

### Create a population DB for all
### data
library(dplyr)
library(tidyr)
library(tidyverse)

steph <- nyando1
steph$count <- 1
steph1 <- steph[, c("IntervieweeName", 
    "Block", "type_r", "area_hectare", 
    "count")]
steph2 <- aggregate(. ~ IntervieweeName + 
    Block + type_r, steph1, sum)

# Create 5 DBH Classes of five DBH
# classes of [<10, 10-20, 20-30,
# 30-40, and >40 cm] <10cm
nyando1_d1 <- filter(nyando1, DBH < 
    10)
nyando1_d1$count <- 1
nyando_d11 <- nyando1_d1[, c("IntervieweeName", 
    "Block", "type_r", "area_hectare", 
    "count")]
aggdata1 <- aggregate(. ~ IntervieweeName + 
    Block + type_r, nyando_d11, sum)
aggdata1$land <- nyando.socio$land_owned_ha[match(aggdata1$IntervieweeName, 
    nyando.socio$name)]
aggdata1$tree_d <- aggdata1$count/aggdata1$land

## Omit infinite values in the
## analyses
aggdata1.b <- aggdata1 %>% filter(tree_d != 
    Inf)

## Lower vs Middle ANOVA
aov.1 <- aov(tree_d ~ Block, data = aggdata1.b)
summary(aov.1)

oneway.test(tree_d ~ Block, data = aggdata1.b)  #ANOVA Unequal var


## Control vs Project ANOVA (Relaxing
## the homogeneity of variance
## assumption)
oneway.test(tree_d ~ type_r, data = subset(aggdata1.b, 
    Block == "Middle"))
oneway.test(tree_d ~ type_r, data = subset(aggdata1.b, 
    Block == "Lower"))

# CALCULATE TREE
# DENSITIES###############################################

## Tree density by
## region#################
aggdata1.b %>% group_by(Block) %>% summarise(mean(tree_d))

## Tree density by treatment group
## (Project vs Control)####
aggdata1.b %>% group_by(type_r) %>% 
    summarise(tree_density = mean(tree_d))

### Tree density for DBH less than
### 10cm##########
nyando_d1 <- filter(nyando1, DBH < 10)
nyando_d1$count <- 1
nyando_d1.1 <- nyando_d1[, c("IntervieweeName", 
    "Block", "type_r", "hectare", "count")]
aggdata1 <- aggregate(. ~ IntervieweeName + 
    Block + type_r, nyando_d1.1, sum)
aggdata1$ha <- steph2$hectare[match(aggdata1$IntervieweeName, 
    steph2$IntervieweeName)]
aggdata1$tree_d <- aggdata1$count/aggdata2$ha

## Lower vs Middle ANOVA
aov.1 <- aov(tree_d ~ Block, data = aggdata1)
summary(aov.1)
oneway.test(tree_d ~ Block, data = aggdata1)  #ANOVA Unequal var

## Control vs Project ANOVA Control
## vs Project ANOVA (Relaxing the
## homogeneity of variance
## assumption)
oneway.test(tree_d ~ type_r, data = subset(aggdata1, 
    Block == "Middle"))
oneway.test(tree_d ~ type_r, data = subset(aggdata1, 
    Block == "Lower"))


### Tree density for DBH
### 10-20cm##########
nyando_d2 <- filter(nyando1, DBH >= 
    10 & DBH <= 20)
nyando_d2$count <- 1
nyando_d21 <- nyando_d2[, c("IntervieweeName", 
    "Block", "type_r", "hectare", "count")]
aggdata2 <- aggregate(. ~ IntervieweeName + 
    Block + type_r, nyando_d21, sum)
aggdata2$ha <- steph2$hectare[match(aggdata2$IntervieweeName, 
    steph2$IntervieweeName)]
aggdata2$tree_d <- aggdata2$count/aggdata2$ha

## Lower vs Middle ANOVA
aov.2 <- aov(tree_d ~ Block, data = aggdata2)
summary(aov.2)
oneway.test(tree_d ~ Block, data = aggdata2)  #ANOVA Unequal var

## Control vs Project ANOVA Control
## vs Project ANOVA (Relaxing the
## homogeneity of variance
## assumption)
oneway.test(tree_d ~ type_r, data = subset(aggdata2, 
    Block == "Middle"))
oneway.test(tree_d ~ type_r, data = subset(aggdata2, 
    Block == "Lower"))

# 20-30CM
nyando_d3 <- filter(nyando1, DBH > 20 & 
    DBH <= 30)
nyando_d3$count <- 1

nyando_d31 <- nyando_d3[, c("IntervieweeName", 
    "Block", "type_r", "hectare", "count")]
aggdata3 <- aggregate(. ~ IntervieweeName + 
    Block + type_r, nyando_d31, sum)
aggdata3$ha <- steph2$hectare[match(aggdata3$IntervieweeName, 
    steph2$IntervieweeName)]
aggdata3$tree_d <- aggdata3$count/aggdata3$ha


## Lower vs Middle ANOVA
aov.3 <- aov(tree_d ~ Block, data = aggdata3)
summary(aov.3)
oneway.test(tree_d ~ Block, data = aggdata3)  #ANOVA Unequal var


## Control vs Project ANOVA
aov.3b <- aov(tree_d ~ type_r, data = aggdata3)
summary(aov.3b)

# TukeyHSD
TukeyHSD(aov.3b)

## Control vs Project ANOVA (Relaxing
## the homogeneity of variance
## assumption)
oneway.test(tree_d ~ type_r, data = subset(aggdata3, 
    Block == "Middle"))
oneway.test(tree_d ~ type_r, data = subset(aggdata3, 
    Block == "Lower"))

# 30-40cm
nyando_d4 <- filter(nyando1, DBH > 30 & 
    DBH <= 40)
nyando_d4$count <- 1
nyando_d41 <- nyando_d4[, c("IntervieweeName", 
    "Block", "type_r", "hectare", "count")]
aggdata4 <- aggregate(. ~ IntervieweeName + 
    Block + type_r, nyando_d41, sum)
aggdata4$ha <- steph2$hectare[match(aggdata4$IntervieweeName, 
    steph2$IntervieweeName)]

aggdata4$tree_d <- aggdata4$count/aggdata4$ha

## Lower vs Middle ANOVA
aov.4 <- aov(tree_d ~ Block, data = aggdata4)
summary(aov.4)
oneway.test(tree_d ~ Block, data = aggdata4)  #ANOVA Unequal var


## Control vs Project ANOVA
aov.4b <- aov(tree_d ~ type_r, data = aggdata4)
summary(aov.4b)

# TukeyHSD
TukeyHSD(aov.4b)

## Control vs Project ANOVA (Relaxing
## the homogeneity of variance
## assumption)
oneway.test(tree_d ~ type_r, data = subset(aggdata4, 
    Block == "Middle"))
oneway.test(tree_d ~ type_r, data = subset(aggdata4, 
    Block == "Lower"))

# >40cm
nyando_d5 <- filter(nyando1, DBH > 40)
nyando_d5$count <- 1

nyando_d51 <- nyando_d5[, c("IntervieweeName", 
    "Block", "type_r", "hectare", "count")]

aggdata5 <- aggregate(. ~ IntervieweeName + 
    Block + type_r, nyando_d51, sum)

aggdata5$ha <- steph2$hectare[match(aggdata5$IntervieweeName, 
    steph2$IntervieweeName)]

aggdata5$tree_d <- aggdata5$count/aggdata5$ha

## Lower vs Middle ANOVA
aov.5 <- aov(tree_d ~ Block, data = aggdata5)
summary(aov.5)
oneway.test(tree_d ~ Block, data = aggdata5)  #ANOVA Unequal var

## Control vs Project ANOVA
aov.5b <- aov(tree_d ~ type_r, data = aggdata5)
summary(aov.5b)

# TukeyHSD
TukeyHSD(aov.5b)

## Control vs Project ANOVA (Relaxing
## the homogeneity of variance
## assumption)
oneway.test(tree_d ~ type_r, data = subset(aggdata5, 
    Block == "Middle"))

# One way ANOVA DBH1 Lower vs Middle
db1 <- aov(DBH ~ type_r, data = nyando1_d1)
summary(db1)

## Combined sample of all DBH
nyando_c <- nyando1
nyando_c$count <- 1
nyando_c1 <- nyando_c[, c("IntervieweeName", 
    "Block", "type_r", "hectare", "count")]
aggdatac <- aggregate(. ~ IntervieweeName + 
    Block + type_r, nyando_c1, sum)
aggdatac$ha <- steph2$hectare[match(aggdatac$IntervieweeName, 
    steph2$IntervieweeName)]
aggdatac$tree_d <- aggdatac$count/aggdatac$ha

## Lower vs Middle ANOVA
aov.c <- aov(tree_d ~ Block, data = aggdatac)
summary(aov.c)
oneway.test(tree_d ~ Block, data = aggdatac)  #ANOVA Unequal var

## Control vs Project ANOVA
aov.cb <- aov(tree_d ~ type_r, data = aggdatac)
summary(aov.cb)

# TukeyHSD
TukeyHSD(aov.cb)

# Levene's Test for Homogeneity of
# Variance (center = median)
oneway.test(tree_d ~ type_r, data = aggdatac)



#### AGB Calculations Combined sample
#### of all DBH
nyando_gb <- nyando1
nyando_gb$num <- rep(1:length(nyando1$IntervieweeName))
nyando_agb <- nyando_gb[, c("Block", 
    "type_r", "Agroforestry Practice", 
    "AGB (Mg)", "num", "mg_ha")]
aggdata_gb <- aggregate(. ~ num + Block + 
    type_r + mg_ha + `Agroforestry Practice`, 
    nyando_agb, sum)

aggdata_gb$agb = aggdata_gb$mg_ha
aggdata_gb$agfp = aggdata_gb$`Agroforestry Practice`
aggdata_gb$Block = as.factor(aggdata_gb$Block)
aggdata_gb$type_r = as.factor(aggdata_gb$type_r)

## Lower vs Middle ANOVA for AGB
oneway.test(agb ~ Block, data = aggdata_gb)  #ANOVA Unequal var

## Lower Project AGB
aggdata_gb1 <- filter(aggdata_gb, type_r == 
    "Lower Project")

# Total AGB under Boundary Planting
sum(filter(aggdata_gb1, agfp == "Boundary planting")$mg_ha)

# Total AGB under Riparian buffers
sum(filter(aggdata_gb1, agfp == "Riparian buffers")$agb)


## Lower Control AGB
aggdata_gb2 <- filter(aggdata_gb, type_r == 
    "Lower Control")

## Total AGB under Boundary Planting
sum(filter(aggdata_gb1, agfp == "Boundary planting")$mg_ha)

# Total AGB under Riparian buffers
sum(filter(aggdata_gb1, agfp == "Riparian buffers")$agb)


## Middle Project AGB
aggdata_gb3 <- filter(aggdata_gb, type_r == 
    "Middle Project")


# Middle Control AGB
aggdata_gb4 <- filter(aggdata_gb, type_r == 
    "Middle Control")

### Tables--Table 5 Calculating Mean
### DBH Middle Control
kab = filter(nyando1, type_r == "Middle Control" & 
    agr == "Boundary planting")
library(summarytools)
dfSummary(kab$DBH, round.digits = 3)
### 
kab1 = filter(nyando1, type_r == "Middle Control" & 
    agr == "Woodlot")
dfSummary(kab1$DBH, round.digits = 3)
### 
kab2 = filter(nyando1, type_r == "Middle Control" & 
    agr == "Riparian buffers")
dfSummary(kab2$DBH, round.digits = 3)

### 
kab3 = filter(nyando1, type_r == "Middle Control" & 
    agr == "MPTs")
dfSummary(kab3$DBH, round.digits = 3)

## Middle Project
kab4 = filter(nyando1, type_r == "Middle Project" & 
    agr == "Boundary planting")
dfSummary(kab4$DBH, round.digits = 3)

### 
kab5 = filter(nyando1, type_r == "Middle Project" & 
    agr == "Woodlot")
dfSummary(kab5$DBH, round.digits = 3)

### 
kab6 = filter(nyando1, type_r == "Middle Project" & 
    agr == "MPTs")
dfSummary(kab6$DBH, round.digits = 3)

### 
kab7 = filter(nyando1, type_r == "Middle Project" & 
    agr == "Hedgerows")
dfSummary(kab7$DBH, round.digits = 3)

### Lower Control
kab8 = filter(nyando1, type_r == "Lower Control" & 
    agr == "Boundary planting")
dfSummary(kab8$DBH, round.digits = 3)

### 
kab9 = filter(nyando1, type_r == "Lower Control" & 
    agr == "Woodlot")
dfSummary(kab9$DBH, round.digits = 3)

## 
kab10 = filter(nyando1, type_r == "Lower Control" & 
    agr == "MPTs")
dfSummary(kab10$DBH, round.digits = 3)

### Lower Project
kab11 = filter(nyando1, type_r == "Lower Project" & 
    agr == "Boundary planting")
dfSummary(kab11$DBH, round.digits = 3)

## 
kab12 = filter(nyando1, type_r == "Lower Project" & 
    agr == "Woodlot")
dfSummary(kab12$DBH, round.digits = 3)

## 
kab13 = filter(nyando1, type_r == "Lower Project" & 
    agr == "Riparian buffers")
dfSummary(kab13$DBH, round.digits = 3)

## 
kab14 = filter(nyando1, type_r == "Lower Project" & 
    agr == "MPTs")
dfSummary(kab14$DBH, round.digits = 3)

## 
kab15 = filter(nyando1, type_r == "Lower Project" & 
    agr == "Hedgerows")
dfSummary(kab15$DBH, round.digits = 3)

tree_d <- rep(c("Acacia seyal", "Balanites aegyptiaca", 
                "Delonix regia", "Eucalyptus sp.", 
                "Eurphobia tirucalli", "Grevillea robusta", 
                "Leucaena trichandra", "Senegalia polyacantha", 
                "Senna siamea", "Thevetia peruviana"), 
              each = 2)

prop_d <- as.numeric(c(0, 9, 32, 13, 
                       0, 3, 0, 6, 0, 23, 0, 4, 0, 17, 
                       10, 0, 0, 4, 13, 5))
site_d <- rep(c("Lower Project", "Lower Control"), 
              10)
## Create dataframe

fig3 <- data.frame(cbind(tree_d, prop_d, 
                         site_d), row.names = NULL)
colnames(fig3) <- c("species", "proportion", 
                    "site")


## Middle Nyando
tree_b <- rep(c("Acacia abyssinica", 
                "Acacia senegal", "Balamites aegyptiaca", 
                "Blighia unijugata", "Carica papaya", 
                "Casuarina equisetifolia", "Croton macrostachyus", 
                "Croton megalocarpus", "Cupressus lusitanica", 
                "Erythrina abyssinica", "Eucalyptus sp.", 
                "Euclea divinorum", "Euphorbia tirucalli", 
                "Ficus sur", "Gliricidia sepum", 
                "Jacaranda mimosifolia", "Kigelia africana", 
                "Mangifera indica", "Markamia lutea", 
                "Persea americana", "Senegalia polyacantha", 
                "Senna siamea", "Spathodea campanulata", 
                "Terminalia brownii", "Vangueria apiculata"), 
              each = 2)

prop_b <- as.numeric(c(10, 5, 0, 2, 
                       5, 0, 0, 1, 2, 0, 0, 0.5, 3, 14, 
                       0, 1, 4, 8, 0, 0.5, 10, 6, 4, 0, 
                       0, 2, 8, 9, 0, 2, 6, 4, 4, 5, 3, 
                       1, 4.5, 4, 1.5, 2, 6, 4, 8, 0, 0, 
                       1, 2.5, 3, 0, 3))

site_b <- rep(c("Middle Project", "Middle Control"), 
              25)

## Create data frame
fig3b <- data.frame(cbind(tree_b, prop_b, 
                          site_b), row.names = NULL)

colnames(fig3b) <- c("species", "proportion", 
                     "site")

## Combine two datasets into one for
## Nyando Add region to Lower
fig3$region <- "Lower"
# Add region to Middle
fig3b$region <- "Middle"

# Combine the two frames
fig3_full <- rbind(fig3, fig3b)
# make proportion numeric
fig3_full$proportion = as.character(fig3_full$proportion)
fig3_full$proportion = as.numeric(fig3_full$proportion)

## Order specied for plotting
facta <- c("Acacia abyssinica", "Acacia senegal", 
           "Balamites aegyptiaca", "Blighia unijugata", 
           "Carica papaya", "Casuarina equisetifolia", 
           "Croton macrostachyus", "Croton megalocarpus", 
           "Cupressus lusitanica", "Erythrina abyssinica", 
           "Eucalyptus sp.", "Euclea divinorum", 
           "Euphorbia tirucalli", "Ficus sur", 
           "Gliricidia sepum", "Jacaranda mimosifolia", 
           "Kigelia africana", "Mangifera indica", 
           "Markamia lutea", "Persea americana", 
           "Senegalia polyacantha", "Senna siamea", 
           "Spathodea campanulata", "Terminalia brownii", 
           "Vangueria apiculata", "Acacia seyal", 
           "Balanites aegyptiaca", "Delonix regia", 
           "Eucalyptus sp.", "Eurphobia tirucalli", 
           "Grevillea robusta", "Leucaena trichandra", 
           "Senna siamea", "Thevetia peruviana")
facta1 <- unique(facta)
facta1 <- sort(facta1, decreasing = TRUE)

fig3_full$species = factor(fig3_full$species, 
                           facta1)
fig3_full$proportion = as.numeric(fig3_full$proportion)

## Create Figure 3 Using ggplot
## barplots Creating figure 3 a and b
## Lower 3a
fig3_full_lower = fig3_full %>% filter(region == 
                                           "Lower")
# Use position=position_dodge()#both
b <- ggplot(data = fig3_full_lower, 
            aes(x = species, y = proportion, 
                fill = site)) + geom_bar(stat = "identity", 
                                         position = position_dodge()) + coord_flip() + 
    scale_fill_manual(values = c("#e34a33", 
                                 "#fdbb84"), labels = c(c("Lower Project-TLP", 
                                                          "Lower Control-TLC")))
b

b.1 <- b + labs(x = "Tree and Shrub Species", 
                y = "Proportions (%)", title = "a). Lower") + 
    theme_bw() + theme(axis.title.x = element_text(size = 14, 
                                                   face = "bold"), axis.title.y = element_text(size = 14, 
                                                                                               face = "bold"), legend.position = "right", 
                       axis.text.x = element_text(size = 12))
b.1

# Middle 3b
fig3_full_middle = fig3_full %>% filter(region == 
                                            "Middle")
# Use position=position_dodge()#both
c <- ggplot(data = fig3_full_middle, 
            aes(x = species, y = proportion, 
                fill = site)) + geom_bar(stat = "identity", 
                                         position = position_dodge()) + coord_flip() + 
    scale_fill_manual(values = c("#a8ddb5", 
                                 "#43a2ca"), labels = c(c("Middle Project-TMP", 
                                                          "Middle Control-TMC"))) + theme_bw()
c  #without labels

c.1 <- c + labs(x = "Tree and Shrub Species", 
                y = "Proportions (%)", title = "b). Middle") + 
    theme(axis.title.x = element_text(size = 14, 
                                      face = "bold"), axis.title.y = element_text(size = 14, 
                                                                                  face = "bold"), legend.position = "right", 
          axis.text.x = element_text(size = 12))

c.1  #with labels



#a8ddb5", "#43a2ca"
